library(dplyr)
library(knitr)
library(maptools)
library(TreeSegmentation)
library(ggplot2)
library(rgl)
library(clue)
library(lidR)
knit_hooks$set(webgl = hook_webgl)
opts_chunk$set(warning=F,message=F)

#set color ramp for treeID
col = pastel.colors(200)

1 Load in ground-truth

shps<-list.files("/Users/ben/Dropbox/Weecology/ECODSEdataset/Task1/ITC/",pattern=".shp",full.names = T)
itcs<-lapply(shps,readShapePoly)

names(itcs)<-sapply(itcs,function(x){
  id<-unique(x$Plot_ID)
  })

itcs<-lapply(itcs,function(x){
  proj4string(x)<-CRS("+init=epsg:32617")
  return(x)
})

2 Example Pipeline

2.1 Read in Data

ground_truth<-itcs[[20]]
fname<-get_tile_filname(ground_truth)
tile<-readLAS(paste("/Users/ben/Dropbox/Weecology/NEON/cropped_",fname,sep=""))
tile@crs<-CRS("+init=epsg:32617")
plot(tile)

You must enable Javascript to view this page properly.

2.2 Compute Segmentation

silva<-silva2016(tile=tile,extra=T)
## [1] "Computing Ground Model"
## [1] "Computing Canopy Model"
## [1] "Clustering Trees"
##    user  system elapsed 
##  14.210   0.373  15.379 
## [1] "Creating tree polygons"
dalponte<-dalponte2016(tile=tile,extra=T)
## [1] "Computing Ground Model"
## [1] "Computing Canopy Model"
## [1] "Clustering Trees"
##    user  system elapsed 
##  12.023   0.277  12.992 
## [1] "Creating tree polygons"
li<-li2012(tile=tile,extra=T)
## [1] "Computing Ground Model"
## [1] "Computing Canopy Model"
## [1] "Clustering Trees"
##    user  system elapsed 
##   2.238   0.086   2.411 
## [1] "Creating tree polygons"

2.3 View 3d segmentation

plot(silva$tile,color="treeID",col=col)

You must enable Javascript to view this page properly.

2.3.1 Overlay ground truth and predictions

plot(ground_truth,col='red')
plot(silva$convex,add=T)

2.3.2 CHM versus predicted polygons

chm=canopy_model(silva$tile)
plot(chm,ext=extent(ground_truth))
plot(ground_truth,add=T,col='red')
plot(silva$convex,add=T)

Okay that’s not great, but let’s keep going for the moment.

2.3.3 Compare Methods

Silva v Dalponte

plot(silva$convex)
plot(dalponte$convex,add=T,col=rgb(0,0,255,20,maxColorValue=255))

Li versus watershed

plot(li$convex)

2.4 Assign Trees

Each tree is assigned based on the maximum overlap. Pairwise membership is done using a Hungarian Algorithm. See clue::solve_LSAP.

assignment<-assign_trees(ground_truth,prediction=silva$convex)

2.5 Calculate Intersection over union

#loop through assignments and get jaccard statistic for each assignment pair
statdf<-calc_jaccard(assignment=assignment,ground_truth = ground_truth,prediction=silva$convex)
ggplot(statdf) + geom_histogram(aes(IoU)) + labs(x="Intersection over union") + theme_bw()

mean(statdf$IoU)
## [1] 0.09047462
median(statdf$IoU)
## [1] 0.07199482

3 As a wrapper for one tile, multiple methods

results<-evaluate(ground_truth=itcs[[5]],algorithm = c("silva","dalponte","li","watershed"),path_to_tiles="/Users/ben/Dropbox/Weecology/NEON/cropped_")
## [1] "Computing Ground Model"
## [1] "Computing Canopy Model"
## [1] "Clustering Trees"
##    user  system elapsed 
##   3.144   0.050   3.310 
## [1] "Creating tree polygons"
## [1] "Computing Ground Model"
## [1] "Computing Canopy Model"
## [1] "Clustering Trees"
##    user  system elapsed 
##   3.044   0.061   3.460 
## [1] "Creating tree polygons"
## [1] "Computing Ground Model"
## [1] "Computing Canopy Model"
## [1] "Clustering Trees"
##    user  system elapsed 
##   0.178   0.006   0.188 
## [1] "Creating tree polygons"
ggplot(results,aes(x=IoU,fill=Method)) + geom_histogram(position = position_dodge()) + theme_bw()

results %>% group_by(Method) %>% summarize(mean=mean(IoU),median=median(IoU))
## # A tibble: 3 x 3
##   Method    mean median
##   <chr>    <dbl>  <dbl>
## 1 dalponte 0.148  0.116
## 2 li       0.272  0.283
## 3 silva    0.139  0.129

4 Wrapper across all tiles

system.time(silvadf<-evaluate_all(itcs=itcs,algorithm = "silva",path_to_tiles="/Users/ben/Dropbox/Weecology/NEON/cropped_"))
qplot(silvadf$IoU)
mean(silvadf$IoU)